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S 1 Abstract: We study the bound state spectrum and the conditions for entering a supercritical 

03 2 regime in graphene with strong intrinsic and Rashba spin-orbit interactions within the 

S 3 topological insulator phase. Explicit results are provided for a disk-shaped potential well 

^ 4 and for the Coulomb center problem. 
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^SJ e 1. Introduction 

O 7 The electronic properties of graphene monolayers are presently under intense study. Previous works 

i— h s have already revealed many novel and fundamental insights; for reviews, see [ 1 ,2] . Following the seminal 

> 9 work of Kane and Mele [3], it may be possible to engineer a two-dimensional (2D) topological insulator 
(TI) phase [4] in graphene by enhancing the — usually very weak [5-7] — spin-orbit interaction (SOI) 
in graphene. This enhancement could, for instance, be achieved by the deposition of suitable adatoms 

12 [8]. Remarkably, random deposition should already be sufficient to reach the TI phase [9-11] where 

13 the effective "intrinsic" SOI A exceeds (half of) the "Rashba" SOI A. So far, the only 2D TIs realized 

14 experimentally are based on the mercury telluride class. Using graphene as a TI material constitutes 

15 a very attractive option because of the ready availability of high-quality graphene samples [1] and the 
is exciting prospects for stable and robust Tl-based devices [4], see also [12,13] 

17 In this work, we study bound-state solutions and the conditions for supercriticality in a graphene- 

18 based TI. Such questions can arise in the presence of an electrostatically generated potential well 
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19 ("quantum dot") or for a Coulomb center. The latter case can be realized by artificial alignment of 

20 Co trimers [14], or when defects or charged impurities reside in the graphene layer. Without SOI, the 

21 Coulomb impurity problem in graphene has been theoretically studied in depth [15-20]; for reviews, see 

22 [1,2]. Moreover, for A = 0, an additional mass term in the Hamiltonian corresponds to the intrinsic 

23 SOI A (see below), and the massive Coulomb impurity problem in graphene has been analyzed in 

24 [21-26]. However, a finite Rashba SOI A is inevitable in practice and has profound consequences. In 

25 particular, A ^ breaks electron-hole symmetry and modifies the structure of the vacuum. We therefore 

26 address the general case with both A and A finite, but within the TI phase A > A/2, in this paper. 

27 Experimental progress on the observation of Dirac quasiparticles near a Coulomb impurity in graphene 

28 was also reported very recently [14], and we are confident that the topological version with enhanced 

29 SOI can be studied experimentally in the near future. Our work may also be helpful in the understanding 

30 of spin-orbit mediated spin relaxation in graphene [27]. 

31 The atomic collapse problem for Dirac fermions in an attractive Coulomb potential, V(r) = 

32 —hvpa/r, could thereby be realized in topological graphene. Here we use the dimensionless impurity 

33 strength 

Ze 2 Z 

a= -—-2.2- (1) 

KaVp K 

34 where Z is the number of positive charges held by the impurity, k a dielectric constant characterizing 

35 the environment, and Vp ~ 10 6 m/s the Fermi velocity. Without SOI, the Hamiltonian is not self-adjoint 

36 for a > a c = 1/2, and the potential needs short-distance regularization, e.g., by setting V(r < R) — 

37 —Uv F o>/R with short-distance cutoff R of the order of the lattice constant of graphene [1,2]. Including 

38 a finite "mass" A, i.e., the intrinsic SOI, but keeping A = 0, the critical coupling a c is shifted to [24] 



1 7T 2 

2 + \n 2 (0.2lAR/hv F ) 



+ , 2, n ~*~,+ s , ( 2 ) 



39 approaching the value a c = 1/2 for R — > 0. In the supercritical regime a > a c , the lowest bound state 

40 "dives" into the valence band continuum (Dirac sea). It then becomes a resonance with complex energy, 

41 where the imaginary part corresponds to the finite decay rate into the continuum. Below we show that 

42 the Rashba SOI provides an interesting twist to this supercriticality story. 

43 The structure of this article is as follows. In Sec. 2 we introduce the model and summarize its 

44 symmetries. The case of a circular potential well is addressed in Sec. 3 before turning to the Coulomb 

45 center in Sec. 4. Some conclusions are offered in Sec. 5. Note that we do not include a magnetic field 

46 (see, e.g., [28,29]) and thus our model enjoys time-reversal symmetry. Below, we often use units with 

47 Tl — Vp — 1. 

48 2. Model and symmetries 



49 2.1. Kane-Mele model with radially symmetric potential 



so We study the Kane-Mele model for a 2D graphene monolayer with both intrinsic (A) and Rashba (A) 
51 SOI [3] in the presence of a radially symmetric scalar potential V(r). Assuming that V(r) is sufficiently 
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52 smooth to allow for the neglect of inter-valley scattering, the low-energy Hamiltonian near the K point 

53 (r = +1) is given by 

H = Tcr x p x + ayPy + tA<j z s z + ^(ra x s y - a y s x ) + V(r), (3) 

54 with Pauli matrices cr x ,y,z (s X)V:Z ) in sublattice (spin) space [1]. The Hamiltonian near the other valley 

55 (K' point) follows for r = —1 in Eq. (3). We note that a sign change of the Rashba SOI, A — > —A, does 

56 not affect the spectrum due to the relation H(—X) = s z H(X)s z . Without loss of generality, we then put 

57 A > and A > 0. 

58 Using polar coordinates, it is now straightforward to verify (see also [21]) that total angular 

59 momentum, defined as 

J z = -id* + s z /2 + ra z /2, (4) 

60 is conserved and has integer eigenvalues j. For given j, eigenfunctions of H must then be of the form 

\ / e^a tiji _(r) \ 



3,T= 



ii4> b 



e^a l>j>+ {r) 



61 Next we combine the radial functions to (normalized) four-spinors, 

( a tdA r ) \ 

®iA r ) = 



e i(i+i)0 a w ,_(r) 
^ ie«*6 Wi _(r) J 



(5) 



a i,jA r ) 
V b uA r ) J 



(6) 



62 In this representation, the radial Dirac equation for total angular momentum j and valley index r = ± 

63 reads 

(H j>T -E)<f> jiT (r) = 0, (7) 

64 with Hermitian matrix operators (note that A denotes the intrinsic SOI and not the Laplacian) 



Hj,- 



es where we use the notation 



(A + V 


v (+) 





\ 


vH 


-A + V 


-A 








-A 


-A + V 


v (+) 


V o 





v<-> 


A + V J 


f -A + V 


v ( -i 





-A \ 


v (+) 


A + V 














A + V 




V "A 





v (+) 


-A + V J 



(8) 



r(±) _ J 



r dr 



(9) 



66 One easily checks that Eq. (8) satisfies the parity symmetry relation 



H 



■3>t 



S y II j -j- ^ 



(10) 
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67 Note that this "parity" operation for the radial Hamiltonian is non-standard in the sense that the valley 

68 is not changed by the transformation a y s y , spin and sublattice are flipped simultaneously, and only the 

69 y-coordinate is reversed. (We will nonetheless refer to a y s y as parity transformation below.) A second 

70 symmetry relation connects both valleys, 

Hj- T = a x Hj )T a x . (11) 

71 Using Eq. (10), this relation can be traced back to a time-reversal operation. Equations (10) and (11) 

72 suggest that eigenenergies typically are four-fold degenerate. 

73 When projected to the subspace of fixed (integer) total angular momentum j, the current density 

74 operator has angular component = o x and radial component J r = —ra y for arbitrary j. When 

75 real- valued entries can be chosen in $^ r (r), the radial current density thus vanishes separately in each 

76 valley. We define the (angular) spin current density as = s z a x . Remarkably, the transformation 

77 defined in Eq. (11) conserves both (total and spin) angular currents, while the transformation in Eq. (10) 

78 reverses the total current but conserves the spin current. Therefore, at any energy, eigenstates supporting 

79 spin-filtered counterpropagating currents are possible. However, in contrast to the edge states found in a 
so ribbon geometry [3], these spin-filtered states do not necessarily have a topological origin. 

81 We focus on one K point (r = +) and omit the r-index henceforth; the degenerate r = — Kramers 

82 partner easily follows using Eq. (11). In addition, using the symmetry (10), it is sufficient to study the 

83 model for fixed total angular momentum j > 0. 



84 2.2. Zero total angular momentum 

85 For arbitrary V(r), we now show that a drastic simplification is possible for total angular momentum 
se j = 0, which can even allow for an exact solution. Although the lowest-lying bound states for the 

87 potentials in Sees. 3 and 4 are found in the j = 1 sector, exact statements about what happens for j = 

88 are valuable and can be explored along the route sketched here. 

89 The reason why j = is special can be seen from the parity symmetry relation in Eq. (10). The parity 

90 transformation a y s y connects the ±j sectors, but represents a discrete symmetry of the j = radial 

91 Hamiltonian Hj =0tT [see Eq. (8)] acting on the four-spinors in Eq. (6). Therefore, the j = subspace 

92 can be decomposed into two orthogonal subspaces corresponding to the two distinct eigenvalues of the 

93 Hermitian operator a y s y . This operator is diagonalized by the matrix 



U 



1 

V2 



( 1 -1 \ 

11 

10 1 

V o i -l o j 



(12) 



94 such that 



UcJySyU 



-1 



diag(l, 1,-1,-1). 



(13) 
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95 
96 



100 
101 



102 
103 



\ 



(14) 



In fact, using this transformation matrix to carry out a similarity transformation, Hj i+ = U Hj i+ U 1 , we 
obtain 

f A + V d r j/r \ 

■l/r-d r -A + V - X j/r 
j/r A + V d r 

j/r -l/r-d r -A + X + V J 

For j = 0, the upper and lower 2x2 blocks decouple. Each block has the signature ("parity") a — ± 
corresponding to the eigenvalues in Eq. (13), and represents a mixed sublattice-spin state, see Eqs. (6) 
and (12). 

For parity a = ±, the 2x2 block matrix in Eq. (14) is formally identical to an effective A = 
problem with j = 0, fixed s z = a, and the substitutions 



A A + aX/2, E ->■ E + cr\/2. 



(15) 



This implies that for j = and arbitrary V(r), the complete spectral information for the full Kane-Mele 
problem (with A ^ 0) directly follows from the A = solution. 



104 2.3. Solution in region with constant potential 



105 
106 
107 
108 



109 
110 
111 
112 
113 
114 



115 
116 



We start our analysis of the Hamiltonian (3) with the general solution of Eq. (7) for a region of constant 
potential. Here, it suffices to study V(r) = 0, since E and V enter only through the combination E — V 
in Eq. (8). In Sec. 3, we will use this solution to solve the case of a step potential. 

The general solution to Eq. (7) follows from the Ansatz 



$,(r) 



( dB^i^/pr) \ 

CsBji^/pr) 
V c 4 5 j+1 ( v ^r) J 



(16) 



(i) 



where the q are real coefficients, Bj is one of the cylinder (Bessel) functions, Bj = Jj or Bj = H 
and p denotes a real spectral parameter. In particular, ^Jp is a generalized radial wavenumber. We here 
assume true bound-state solutions with real-valued energy. However, for quasi-stationary resonance 
states with complex energy, p and the c, may be complex as well. 

Using the Bessel function recurrence relation, V ^ Bj (y/pr) = y/pBj T i(^r), the set of four coupled 
differential equations (7) simplifies to a set of algebraic equations, 



( A-E 










\ 


( Cl ) 


Vp 


-A-E 


-A 







c 2 





-A 


-A-E 


Vp 




c 3 


V o 







A-E 


) 


V C 4 ) 



0. 



(17) 



Notably, j does not appear here, and therefore the spectral parameter p depends only on the energy E. 
The condition of vanishing determinant then yields a quadratic equation for p, with the two solutions 



p± = (E — A)(E — E±), E± = -A±\. 



(18) 
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117 
118 
119 
120 



121 
122 



Which Bessel function is chosen in Eq. (16) now depends on the sign of p± and on the imposed regularity 
conditions for r — > and/or r — > oo. 

For p± > 0, a solution regular at the origin is obtained by putting Bj = Jj, which describes standing 
radial waves. Equation (17) then yields the unnormalized spinor 

fp±r) \ 

TJj(y/P~±r) 

/PI 



( v 7 ^ J i 



j,P±>0 



(19) 



For p± < 0, instead it is convenient to set Bj = Hj in Eq. (16). Using the identity H 
—.e~ % ^l 2 Kj^z), the unnormalized spinor resulting from Eq. (17) then takes the form 



( 



-p± 



(20) 



123 
124 



E-A 

Ti^(V=pIr) 

where the modified Bessel function Kj(^—p±r) describes evanescent modes, exponentially decaying at 
infinity. 



125 2.4. Solution without potential 



126 
127 
128 
129 
130 
131 

132 
133 
134 
135 
136 
137 



In a free system, i.e., when V(r) = for all r, the only acceptable solution corresponding to a physical 
state is obtained when p± > [30]. For A < A/2, at least one p± > in Eq. (18) for all E, and the 
system is gapless. However, the TI phase defined by A > A/2 has a gap as we show now. 

For A > A/2, Eq. (18) tells us that for E > A and for E < E_, both solutions p± are positive and 
hence (for given j and r) there are two eigenstates $ JiP± for given energy E. However, within the energy 
window [with E± in Eq. (18)] 

E_<E < E + , (21) 

we have p + > and p_ < 0, i.e., only the eigenstate &j, p+ represents a physical solution. Both p± 
are negative when E + < E < A, and no physical state exists at all. This precisely corresponds to 
the topological gap in the TI phase [3]. Note that due to the Rashba SOI, the valence band edge is 
characterized by the two energies E±, with halved density of states in the energy window (21). One 
may then ask at which energy (E + or EJ) the supercritical diving of a bound state level in an impurity 
potential takes place. 



3. Circular potential well 



139 3.1. Bound states 

ho In this section, we study a circular potential well with radius R and depth V > 0, 

V{r) = { ~ Vo > r<R (22) 
* 0, r > R 
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Figure 1. Bound-state spectrum (Eb) vs Rashba SOI (A) for a circular potential well with 
depth V = 2A and radius R = 3/A. Only the lowest-energy states with j = 0, 1, 2 are 
shown. The red dotted line indicates E + = —A + A. The left panel shows j = bound states 
with parity a = ±. The right panel shows j = 1,2 bound states. The inset displays the j = I 
bound-state energies vs potential depth V for A = 0.6A. At some threshold value Vo = V t 
(where Vt — for the lowest state shown), a new bound state emerges from the conduction 
band. This state dives into the valence band for some critical value V = V c > V t , where the 
valence band edge is at energy E B = E + = — 0.4A. For the second bound state in the inset, 
V t (V c ) is shown as red (blue) triangle. 




L i I i I i I i I L i I i I i l_ i I 

0.5 1 1.5 0.5 1 1.5 2 
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141 We always stay within the TI phase A > A/2, where bound states are expected for energies E = E B 

142 in the window max(A — Vo, E + ) < Eb < A. For r < R, the corresponding radial eigenspinor [see 

143 Eq. (6)] is written with arbitrary prefactors A± in the form 

$<(r) = £A<$^ ± (r), (23) 
± 

144 with Eq. (19) for $j,p ± (r). Here, the p± > follow from Eq. (18) by including the potential shift, 

p± = (E + V - A)(E + V Q - E ± ). (24) 

145 For r > R, the general solution is again written as 

$>(r)=5>>$ JiP± (r). (25) 
± 

146 However, now $j, p ± is given by Eq. (20), since p± < for true bound states with only evanescent states 

147 outside the potential well. 

us The continuity condition for the four-spinor at the potential step, $f(R) = &f(R), then yields a 

149 homogeneous linear system of equations for the four parameters (A± > ). A nontrivial solution is only 

150 possible when the determinant of the corresponding 4x4 matrix C(E) (which is too lengthy to be given 

151 here but follows directly from the above expressions) vanishes, 

det[C(E)\ = 0. (26) 

152 Solving the energy quantization condition (26) then yields the discrete bound-state spectrum (Eb). It is 

153 then straightforward to determine the corresponding spinor wavef unctions. 

154 Numerical solution of Eq. (26) yields the bound-state spectrum shown in Fig. 1. When Vo exceeds a 

155 (j-dependent) "threshold" value, V t , a bound state splits off the conduction band edge. When increasing 

156 Vo further, this bound-state energy level moves down almost linearly, cf. inset of Fig. 1, and finally 

157 reaches the valence band edge E + = —A + A at some "critical" value Vo = V c . (For j = 0, we will 

158 see below that this definition needs some revision.) Increasing V even further, the bound state is then 

159 expected to dive into the valence band and become a finite-width supercritical resonance, i.e., the energy 

160 would then acquire an imaginary part. 

161 3.2. Zero angular momentum states 

162 Surprisingly, for j = 0, we find a different scenario where supercritical diving, with finite lifetime of 

163 the resonance, happens only for half of the bound states entering the energy window (21). Noting that 

164 states with different parity a = ± do not mix, see Sec. 2.2, we observe that all a = + bound states enter 

165 the valence band as true bound states (no imaginary part) throughout the energy window (21) while the 

166 valence band continuum is spanned by the a = — states. We then define V c for (j = 0, a = +) bound 

167 states as the true supercritical threshold where E B = E_ = —A — A. However, the (j = 0, a = —) 

168 bound states become supercritical already when reaching E + = — A + A. 

169 Therefore an intriguing physical situation arises for j = in the energy window (21). While a = + 

170 states are true bound states (no lifetime broadening), they coexist with a = — states which span the 
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171 valence band continuum or possibly form supercritical resonances. For E < E_, however, all bound 

172 states dive, become finite-width resonances, and eventually become dissolved in the continuum. 

173 3.3. Threshold for bound states 

174 Returning to arbitrary total angular momentum j, we observe that whenever V hits a possible 

175 threshold value V t , a new bound state is generated, which then dives into the valence band at another 

176 potential depth V = V c (and so on). Analytical results for all possible threshold values V t follow by 

177 expanding Eq. (26) for weak dimensionless binding energy 6=1 — EbJ A. For 5 <C 1 and j = 1, 

178 Eq. (26) yields after some algebra 

A = WVF/R) 2 e-* ( I - A V /2 e -^ E± , ±Jo fa) Mfe ) ) (2?) 



A 2 V1-A 2 VI + A 
z ± = Jv Q (2A ± A + V Q )R, 



179 where 7 w 0.577 is the Euler constant and A = A/2A. The binding energy approaches zero for V$ — > 0, 

180 where Eq. (27) simplifies to 



2(hv F /R) 2 (1~X\ 



r \ A/2 „ 

e 7 v ^ . (28) 



A 2 \/l-A 2 



181 For vanishing Rashba SOI A = 0, this reproduces known results [25]. For any A < 2 A, we observe that 

182 the j — 1 bound state in Eq. (28) exists for arbitrarily shallow potential depth Vo. 

183 The threshold values V t for higher-lying j = 1 bound states also follow from the binding energy (27), 

184 since 6 vanishes for J\{z + ) = and for J\(z_) = 0. When one of these two conditions is fulfilled 

185 at some Vo = V t , a new bound state appears for potential depth above V t . This statement is in fact 

186 quite general: By similar reasoning, we find that the threshold values V t for j = follow by counting the 

187 zeroes of Jo (z± ) . Without SOI, this has also been discussed in [3 1 ] . Note that this argument immediately 

188 implies that no bound state with j = exists for Vo — > 0. 

189 From the above equations, we can then infer the threshold values V t for all bound states with j = 

190 or j = 1 in analytical form. These are labeled by n = 1,2,... and a = ± (for j = 0, a corresponds to 

191 parity), 

V tj> ,± = (A ± A/2) [-1 + ^l + 1 lJ[R(A±X/2)} 2 ] , (29) 

192 where 7^ is the nth zero of the Jj Bessel function. 

193 Likewise, for j > 1, the condition for the appearance of a new bound state is 

J,{z T ) = 0. (30) 



E 

± 



194 Close examination of this condition shows that no bound states with j > 1 exist for V — > 0. We 

195 conclude that bound states in a very weak potential well exist only for j = 1. 
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Figure 2. Critical potential depth V c for the lowest j = I bound state level in a disk with 
RA = 3. The obtained A = value matches the analytical prediction V c ~ 2.28A from 
Eq. (31), while V c — > near the border of the TI phase (A — > 2 A). Inset: V c vs radius R with 
several values of A (given in units of A) for the lowest bound state. 



2 
1.5 

< 

^ i 

0.5 



w 0.5 1 1.5 2 

A/A 

196 3.4. Supercritical behavior 

197 As can be seen in Fig. 1, the lowest j = 1 bound state is also the first to enter the valence band 

198 continuum for V$ = V c . For A = 0, the critical value is known to be [25] 

K = A(l + v /l + 7o 2 1 /[i?A]2). (31) 

199 with 7 0) i ~ 2.41. The energy of the resonant state acquires an imaginary part for Vq > V c [25]. For 

200 A > 0, we have obtained implicit expressions for V c , plotted in Fig. 2. Note that these results reproduce 

201 Eq. (31) for A — > 0. The almost linear decrease of V c with increasing A, see Fig. 2, can be rationalized by 

202 noting that the valence band edge is located at E + = — A + A. Thereby supercritical resonances could be 

203 reached already for lower potential depth by increasing the Rashba SOI. Similarly, with increasing disk 

204 radius R, the critical value V c decreases, see the inset of Fig. 2. For the lowest (j = 0, o = ±) bound 

205 state, the critical value in fact follows in analytical form, 

V c>a=± = (A ± A/2) (l + ^l + 7l 2 1 /[(A±A/2) J R]^ , (32) 

206 where 711 w 3.83. 

207 Since the parity decoupling in Sec. 2.2 only holds for j = 0, it is natural to expect that all j 7^ bound 

208 states turn into finite-width resonances when Eb < E + . This expectation is confirmed by an explicit 

209 calculation as follows. Within in the window E_ < E B < E + , a true bound state should not receive 

210 a contribution from $ JiP+>0 (r) for r > R, but instead has to be obtained by matching an Ansatz as in 
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213 



Figure 3. Bound state energies with angular momentum j = (E B in units of A) vs 
dimensionless impurity strength a for the Coulomb problem with regularization parameter 
RA = 0.01 and Rashba SOI A = 0.6A. Solid black (dashed blue) curves correspond to 
parity a = + (a = — ). Results for radial number n = 1, 2, 3 (with increasing energy) are 
shown. Red dotted lines denote E = E±. 




Eq. (23) for the spinor state inside the disk (r < R) to an evanescent spinor state oc $j )P _<o(r > R). 
However, the matching condition is then found to have no real solution E B , i.e., there are no true bound 
states with j '> =/ in the energy window (21). We therefore conclude that all j / bound states turn 

214 supercritical when Eb < E + . Note that this statement includes the lowest-lying bound state (which has 

215 j = 1). This implies that a finite Rashba SOI can considerably lower the potential depth V c required for 

216 entering the supercritical regime. 



4. Coulomb center 



220 



We now turn to the Coulomb potential, V(r) = —a/r, generated by a positively charged impurity 
located at the origin, with the dimensionless coupling strength a in Eq. (1). We consider only the TI 
phase A > A/2 and analyze the bound-state spectrum and conditions for supercriticality. Again, without 
loss of generality, we focus on the K point only (r = +), and first summarize the known solution for 
A = [2,21,24]. In that case, s z = ± is conserved, and the spin-degenerate bound-state energies are 
labeled by the integer angular momentum j and a radial quantum number n = 1, 2, 3, . . . (for j > 0, 
n = is also possible), 

-1/2 



E hn {\ = 0) = A 1 



a 



(n + - 1/2)2 _ a 2j 



(33) 



The corresponding eigenstates then follow in terms of hypergeometric functions. The lowest bound state 
is Ej=i jn=0 = Ay/1 — 4a 2 , which dives when a = a c = 1/2; note that a c precisely corresponds to V c in 
Sec. 3. In particular, for (j — 0, a — ±) states we define a c in the same manner. Next we discuss how 



228 this picture is modified when the Rashba coupling A is included. 
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230 



231 



232 



236 



237 



238 



243 



Figure 4. Main panel: Critical Coulomb impurity strength a c vs Rashba SOI A for RA = 
0.01 and the lowest (J — 0, a = ±) bound states. Inset: a c vs cutoff scale R for A = 0.6A. 




Following the arguments in Sec. 2.2 for j = 0, the combination of Eq. (33) with Eq. (15) immediately 
yields the exact bound-state energy spectrum (n — 1, 2, 3, . . .), 



/ 



-1/2 



-Ej=0,n,o-=± 



(A ± A/2) 



1 + 



a 



T X/2. 



(34) 



The corresponding eigenstates then also follow from Refs. [21,24]. The very same reasoning also applies 
to a regularized 1/r potential [23,24], where V(r < R) is replaced by the constant value V = —a/R. 
Here, R is a short-distance cutoff scale of the order of the lattice spacing. The solution of the bound-state 
problem then requires a wavefunction matching procedure, which has been carried out in Ref. [24]. 
Thereby we can already infer all bound states for j = 0. 

Figure 3 shows the resulting j = bound-state spectrum vs a for the regularized Coulomb potential. 
Within the energy window (21), we again find that states with parity a = + remain true bound states 
that dive only for Eb < E_, while a = — states show supercritical diving already for Eb < E + . Figure 
4 shows the corresponding critical couplings a c for a = ±, where the lowest j = bound state with 
parity a turns supercritical. Note that for finite R and A — > 0, a unique value for a c is found, while 
for A / two different critical values for a c are found. However, this conclusion holds only for finite 
regularization parameter R, i.e., it is non-universal. As seen in the inset of Fig. 4, in the limit R — > 0, 
both critical values for a c approach a c — 1/2 again, which is the value found without SOI. 

Finally, for j ^ 0, we can then draw the same qualitative conclusions as in Sec. 3.4 for the potential 
well. In particular, we expect that all j ^ bound states turn supercritical when their energy Eb reaches 
the continuum threshold at Eb = E + = — A + A. 



5. Conclusions 
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248 In this work, we have analyzed the bound-state problem for the Kane-Mele model of graphene with 

249 intrinsic (A) and Rashba (A) spin-orbit couplings when a radially symmetric attractive potential V(r) is 

250 present. We have focussed on the most interesting "topological insulator" phase with A > A/2. The 

251 Rashba term A leads to a restructuring of the valence band, with a halving of the density of states in the 

252 window < E < E+, where E± = —A ± A. This has spectacular consequences for total angular 

253 momentum j = 0, where the problem can be decomposed into two independent parity sectors (a = ±). 

254 The a = + states remain true bound states even inside the above window and coexist with the continuum 

255 solutions as well as possible supercritical resonances in the a = — sector. However, all j ^ bound 

256 states exhibit supercritical diving for E < E + , where the critical threshold (V c or a c for the disk or the 

257 Coulomb problem, respectively) is lowered when the Rashba term is present. We hope that these results 

258 will soon be put to an experimental test. 
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